Large well-relaxed models of vitreous silica, 
coordination numbers and entropy 
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A Monte Carlo method is presented for the simulation of vitreous silica. Well-relaxed networks 
of vitreous silica are generated containing up to 300,000 atoms. The resulting networks, quenched 
under the BKS potential, display smaller bond-angle variations and lower defect concentrations, as 
compared to networks generated with molecular dynamics. The total correlation functions T{r) of 
our networks are in excellent agreement with neutron scattering data, provided that thermal effects 
and the maximum inverse wavelength used in the experiment are included in the comparison. A 
procedure commonly used in experiments to obtain coordination numbers from scattering data is 
to fit peaks in rT(r) with a gaussian. We show that this procedure can easily produce incorrect 
results. Finally, we estimate the configurational entropy of vitreous silica. 

PACS numbers; 61.43.Fs,61.43.Bn,61.12.Bt 
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I. INTRODUCTION 

If liquid silica (Si02) is cooled below its melting tem- 
perature, it usually does not crystallize, but stays a su- 
percooled liquid for an extended period of time. Upon 
further cooling, to below a certain temperature known 
as the glass transition temperature Tg, it becomes highly 
viscous and shows many properties of ordinary solids, 
but stays disordered and does not evolve to the thermo- 
dynamically stable crystalline phase, at least on earthly 
time scales. The study of silica in its disordered phase 
(vitreous silica) is hampered by a lack of knowledge on its 
microscopic structure. Experimental techniques alone, 
such as neutron and X-ray scattering, which can uniquely 
resolve the structure of a crystal, are less successful in this 
case because of the lack of a repeating unit cell. Typi- 
cally, only information averaged over many atoms, such 
as the radial distribution function, can be obtained from 
experiments. 

An alternative method to study vitreous silica is 
through computer simulations. The usual computational 
approach to generate vitreous silica structures is to sim- 
ulate a quench from the melt within the framework of 
molecular dynamics (MD) . This resembles closely the ex- 
perimental process used to prepare vitreous silica, except 
that the typical computational cooling rates are about 
ten orders of magnitude faster than experimental ones. 
In agreement with experiment, these MD simulations re- 
sult in structures of vitreous silica in which almost all 
silicon atoms are bonded to four oxygen atoms and al- 
most all oxygen atoms are bonded to two silicon atoms, 
without any sign of long-range order |l) . This in turn is in 
agreement with the continuous random network (CRN) 
model, proposed by Zachariasen H. 



The high cooling rates used in MD simulations result 
in structures with a strain higher than observed in exper- 
iments. This shows up, for instance, as an anomalously 
high density of coordination defects and a larger spread in 
bond lengths and bond angles. In this paper we present 
an alternative approach to generate models of vitreous 
silica. In contrast to MD, we do not attempt to simulate 
the details of the dynamics of the melt-quenching pro- 
cess. Instead, we use a Monte Carlo (MC) scheme with 
an artificial dynamics consisting of bond transpositions. 
Our scheme is similar in spirit to the one used by Tu 
et al. in their study of silica and silica sub-oxides ^, Q, 
which in turn is based on the algorithm of Wooten, Winer 
and Weaire (WWW) for the generation of four-fold co- 
ordinated CRNs |, 1. 

The outline of this paper is as follows. We first de- 
scribe the MC scheme used by Tu et al. We then move 
on to describe a number of optimizations to this algo- 
rithm. These optimizations allow us to generate large 
and well-relaxed silica networks containing up to 300,000 
atoms. The properties of these networks are discussed 
and compared to MD-prepared networks and experiment 
in Section Vl . We end with a number of conclusions in 



Section VIII 



II. 



BOND-SWITCHING ALGORITHM FOR 
SILICA 
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In the approach followed by Tu and co-workers, a vit- 
reous silica network consists of Ng silicon atoms and 
No — 2Ns oxygen atoms. The total number of atoms 
in the network thus equals N = Ns + No- At all times, 
an explicit list of bonds is maintained which uniquely 
determines the connectivity or topology of the network. 
In the list of bonds, each silicon atom is bonded to four 
oxygen atoms, and each oxygen atom is bonded to two 
silicon atoms. The total number of bonds in the list thus 
equals ANs. 
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To calculate energy and forces the following potential 
is used: 

i 

+ i ^ fee (cos cos 6*0)^, (1) 

where the first summation runs over the list of bonds, and 
the second summation over all pairs of bonds that share 
one atom; bi represents the length of the i-th bond, and 
9ij is the angle between bonds i and j (which share one 
atom). The potential parameters for the two-body Si-0 
interaction are given by kb = 27.0 eV and bo — 1-60 A. 
For the three-body 0-Si-O interaction common values are 
kg = 4.32 eV and cos^o = —1/3; for Si-O-Si interactions 
these values are kg = 0.75 eV and cos^o = —1- Since 
the number of bonds in the network grows linearly with 
the number of atoms it contains, each energy evaluation 
is an 0{N) operation. 

The network is evolved by making explicit changes in 
the list of bonds, each time followed by a local energy 
minimization with respect to the geometric coordinates 
of all atoms (using for example the method of steepest 
descent). The changes in the list of bonds consist of 
bond transpositions as illustrated in Fig. |^. Each bond 
transposition is accepted with the Metropolis acceptance 
probability 0: 



P„ = min 



1, exp 



- Ef 

keT 



(2) 



where fc^ is the Boltzmann constant, T is the tempera- 
ture, and Eb and Ef are the (minimized) energies of the 
network before and after the proposed bond transposi- 
tion. 



III. OPTIMIZED BOND-SWITCHING 
ALGORITHM FOR SILICA 



the total energy is harmonic around the minimum with 
curvature c, the energy E{Rq) at the (unknown) mini- 
mum i?o can be estimated from a nearby configuration 
R as E{Ra) ~ E{R) - |^|Vc, with |F| the magnitude 
of the total force at R. Minimization is aborted as soon 
as it becomes clear that the threshold energy will not be 
reached and the bond transposition will be rejected. This 
leads to a large reduction in the number of force evalua- 
tions associated with rejected bond transpositions. 

The second optimization is aimed at exploiting the lo- 
cal nature of the bond transposition depicted in Fig. ^. 
Immediately after a bond transposition, only a small 
cluster of atoms in the network will experience a sig- 
nificant force. This cluster consists of the atoms 
directly involved in the bond transposition, marked 
{S'l, 52, 01, 02, 03} in Fig. |l|, and of nearby atoms (typ- 
ically extending to the fourth neighbor shell around these 
five transposition atoms). The number of atoms in such 
a cluster is about 400. It therefore suffices to calculate 
the force locally (only for the 400 or so atoms inside the 
cluster) rather than globally (for all the atoms in the net- 
work). Calculating the force on a cluster of atoms is an 
0(1) operation, which means that it is independent of the 
total number of atoms in the network. Local force cal- 
culations are therefore much cheaper than global 0{N) 
force calculations. 

Already after a converged local minimization, it is of- 
ten clear that the threshold energy will not be reached, 
and the bond transposition can be rejected. If this is not 
the few global minimization steps are required 

additionally, usually resulting in an accepted bond trans- 
position. The combination of these two improvements 
reduces the computational effort from 0{N) operations 
per attempted bond transposition to 0(1) operations per 
attempted bond transposition, plus a few 0{N) opera- 
tions per accepted bond transposition. Since common 
acceptance probabilities are well below 1%, the speed-up 
for large systems is about two orders of magnitude. 



We have optimized the bond-switching algorithm of Tu 
et al. for the generation of large and well relaxed silica 
networks. These optimizations are similar in spirit to the 
optimizations used in the scalable WWW algorithm ||, 

I- 

The first optimization is aimed at reducing the CPU 
time spent on rejected bond transpositions. After a bond 
transposition in the original algorithm, the energy of the 
network is always completely minimized. After the min- 
imization, the bond transposition is either accepted or 
rejected based on the Metropolis probability. 

In contrast, we reformulate the Metropolis algorithm. 
We first determine a threshold energy Et ~ E — 
kBT\n{l — r), where r is a random number uniformly 
drawn from the interval [0, 1). We then proceed with the 
minimization procedure. During minimization, the con- 
verged energy is continuously estimated. Assuming that 



IV. INITIAL CONFIGURATIONS 

The optimized bond-switching algorithm for silica re- 
quires an initial network, which should be free of co- 
ordination defects and crystalline regions. As noted in 
Ref. |l^, four-fold coordinated CRNs decorated with one 
oxygen atom on each bond already provide structures 
that compare reasonably well to vitreous silica. As ini- 
tial network we therefore use a periodic, four-fold coor- 
dinated CRN generated as described in Refs. ^. This 
CRN serves as the silicon backbone. Next, one oxygen 
atom is placed on the center of each Si-Si bond to obtain 
a properly coordinated silica network. The simulation 
volume is then scaled to obtain the experimental density 
of silica p = 2.20 g cm'^ [|l). 
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FIG. 1: Bond transposition used in the generation of vitreous silica structures. Two silicon atoms {SI, 52} and three oxygen 
atoms {Ol, 02, 03} are selected following the geometry shown left. Next, bonds 5*1 — Ol and S2 ~ 03 are broken and two 
new bonds SI — OS and S2 — Ol are created to obtain the geometry shown right. After the bond transposition, the geometric 
coordinates of the atoms are relaxed (not shown). 



V. BKS SAMPLING 

A computationally more expensive, but also more re- 
alistic potential for silica is the BKS potential ||12i 



Ebks 




Aij exp {-Bijrij) 




(3) 



where the sum runs over all pairs of atoms in the sys- 
tem, Vij is the distance between atoms i and j, and qi is 
the charge of atom i. The atomic charges and the val- 
ues of the potential parameters Aij , Bij and Cij are also 



listed in Ref. 12, We use the fast multipole code devel- 
oped at Duke University to compute the Coulomb 
and van der Waals interactions. The exponential term 
is truncated at 5.50 A and then shifted, as described in 
Ref. I. 

The bond-switching algorithm described above gener- 
ates networks that are local energy minima according to 
the potential of Tu et al. The algorithm can be changed 
such that these networks will be sampled from the Boltz- 
mann distribution according to the BKS potential at 
temperature T. To this end, the acceptance probabil- 
ity of Eq. (0) is replaced by: 



P = min 



X mm 



l,exp 
1, cxp 



AE 



'TT 



AE 



BKS 



V keT' ksT 



(4) 



where AEtt and AEbks are the energy differences be- 
tween the networks before and after the bond transposi- 
tion, calculated with the potential of Tu et al. and the 
BKS potential, respectively. The algorithm satisfies de- 
tailed balance, also when the fictitious temperatures T 
and T' are different. In the expression above, the first 
factor biases moves according to the potential of Tu et 
al. at temperature T', after which the second factor cor- 
rects for the difference between this potential at temper- 
ature T' and the BKS potential at temperature T . 

With the above probability, bond transpositions that 
are accepted with the potential of Tu et al. undergo an 



additional accept/reject decision based on the BKS po- 
tential. By increasing the temperature T' we can increase 
the number of moves that will be accepted under the po- 
tential of Tu et al. In the limit T' oo, all such moves 
are accepted. In this limit, the accept/reject decision is 
based solely on the BKS potential and T. While correct, 
this situation is undesirable because BKS energy evalu- 
ations are expensive. In practice, a lower T' is used and 
the temperatures {T, T'j are tuned to maximize the de- 
crease of the BKS energy of the network as a function of 
simulation time. 



VI. RESULTS 

We have used the optimized bond-switching algorithm 
described in Section |l| to generate vitreous silica net- 
works containing 3000 atoms ('3k'), 60,000 atoms ('60k') 
and 300,000 atoms ('300k'). For the 3000-atom and the 
60,000-atom networks the bond-switching algorithm was 
used in conjunction with BKS sampling at temperatures 
kgT = 0.30 eV and keT' = 0.70 eV (see Section 0). 
For the 300,000-atom network, BKS sampling proved to 
be too computationally demanding and was not imple- 
mented; for this network the Metropolis probability of 
Eq. (^) was used with ksT = 0.15 eV and the potential 
of Tu et al. 

All three networks were evolved with approximately 
10 attempted bond transpositions per atom from their 
starting configurations. The energy of the resulting net- 
works was then minimized in a single quench using the 
method of steepest descent and the BKS potential with- 
out volume optimization. We observed that if volume 
optimization were used with the BKS potential, the den- 
sity of the networks becomes unphysically high, typically 
in the range 2.27-2.37 g cm^"^. Similar high densities are 
also reported in a number of MD studies I], |l|. 

In Fig. H we show the radial distribution func- 
tions (RDFs) gsiSi{r), 9Sio{r), and goo{r) for model 
'60k'. The RDFs are density normalized such that 
limr^oo <7q/3('') = 1- The RDFs for models '3k' and 
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FIG. 2: Species-dependent density-normalized radial distri- 
bution functions gapir) for network '60k' after relaxation at 
zero temperature with the BKS potential p^]. 



'300k' are similar (not shown). 

The Si-0 RDF was used to extract a cutoff for the Si- 
O bond length at Vc = 1-80 A, where Tc is taken at the 
minimum between the first and the second peak. Based 
on this Tc we find structural properties as reported in 
Table ^. For comparison, we also report results obtained 
with MD using the same BKS potential. Compared to 
MD, we observe that our networks are better relaxed, as 
is evident from the smaller variations in the 0-Si-O and 
Si-O-Si bond angles. Moreover, our networks are nearly 
free of coordination defects. 

A more stringent test is to compare the properties of 
our networks to experimental data. The crucial quan- 
tity to consider is the total correlation function T(r) as 



TABLE I: Structural and energetic properties of our vitre- 
ous silica networks at zero temperature, after a local energy 
minimization under the BKS potential at the experimental 
density. For comparison we also report properties of silica 
networks prepared by a quench to zero temperature under 
MD using the same potential. The number of atoms in 
each structure is given by A''; E is the BKS energy per silicon 
atom in eV; p is the density in g cm~^; the Si-0 bond length 
and its rms variation are given in A; the mean bond angles 
and their variations (rms and FWHM) are given in degrees; 
Z4 and Z2 are the percentages of perfectly coordinated silicon 
and oxygen atoms, respectively, based on a Si-O bond-length 
cut-off at 1.80 A. 





'3k' 


'60k' 


'300k' 


MD 


N 


3000 


60,000 


300,000 


1002 


E 


-58.12 


-58.10 


-58.09 




P 


2.20" 


2.20" 


2.20" 


2.27-2.37 


Si-0 










mean 


1.606 


1.608 


1.606 


1.595' 


rms 


0.010 


0.011 


0.011 




O-Si-0 










mean 


109.44 


109.43 


109.43 


108.3 


rms 


3.95 


4.59 


4.32 




FWHM 


8.3 


9.3 


9.8 


12.8 


Si-O-Si 










mean 


153.89 


153.57 


153.00 


152 


rms 


11.75 


11.72 


11.94 




FWHM 


34 


33.3 


34.5 


35.7 


24 


100% 


100% 


99.997% 


99.8% 


22 


100% 


100% 


99.998% 


99.8% 



"Fixed to the experimental density p = 2.20 g 
'Location of the first peak in the Si-O RDF. 
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obtained in X-ray or neutron scattering experiments. In 
these experiments, the quantity that can be measured di- 
rectly is the total interference function Qi{Q); the corre- 
lation function T{r) is obtained via a Fourier transform of 
Qi{Q). In practice, the resolution in T{r) is determined 
by the maximum attainable inverse wavelength Qmax in 
the experiment. As discussed in Ref. |l^, peaks in T(r) 
are significantly broadened because Qmax is not infinite 
in real experiments. In addition to Fourier broadening, 
there is thermal broadening because experiments are typ- 
ically carried out at room temperature. For a meaning- 
ful comparison it is essential to take these factors into 
account. 

To capture the effect of thermal broadening, all three 
networks are thermalized over approximately 1 ps with 
MD using the BKS potential in the NVE-ensemble at 
T K, 300 K. The structural properties of the thermalized 
networks are summarized and compared to experiment 
in Table ||. We observe good agreement on the Si-0 
bond length, the variation in the bond length and the 
mean bond angles. Compared to experiment, our net- 
works slightly overestimate the variation in the 0-Si-O 
bond angle. This may indicate that experimental vitre- 
ous silica is more relaxed than our networks. Experimen- 
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TABLE II; Structural properties of networks '3k', '60k' and 
'300k' after thermalization at 300 K, based on a Si-O bond- 
length cut-off at 1.80 A. Shown are the mean Si-O bond length 
with rms variation in A; the mean bond angles with rms varia- 
tion and FWHM of the corresponding distribution in degrees; 
and the mean coordination number z of silicon atoms. For 
comparison we also report values obtained in experiments. 





'3k' '60k' '300k' Experiment 


Si-O 

mean 

rms 


1.610 1.614 1.612 1.609 O] 
0.042 0.042 0.041 0.047 gj 


O-Si-O 

mean 

rms 

FWHM 


109.38 109.37 109.37 109.7 ± 0.6 |||] 
5.79 6.27 6.31 4.5 @ 
13.1 14.4 14.5 


Si-O-Si 

mean 

rms 

FWHM 


153.01 152.74 152.20 148 ||l|l;152 [|| 
12.22 12.31 12.60 7.5 
35 34.4 36.2 12.8 p|; 26 p|; 36 |T7|] 


z 


4.000 4.000 4.000 



tal estimates for the variation in the Si-O-Si bond angle 
range from 12.8 degrees to 36 de grees [0. Compared 
to these data our networks coincide with the higher value. 
The coordination numbers reported in Table O indicate 
that the vast majority of atoms is still properly coordi- 
nated after thermalization. 

Next, Fourier broadening is taken into account. In 
the absence of Fourier broadening, the total correlation 
function of a network is simply a weighted sum of the 
RDFs: 



Tir) 



{wsiSi X gsiSi{r) 



WSiSi + WSiO + Woo 

+ wsio X gsto{r) + woo x goo{r)] , 



(5) 



in which each RDF is normalized such that 
limr^oo T{r) — r. We focus on neutron scattering, 
in which case the weights are given by: wsiSi = 0.1722, 
wsio = 2 X 0.4817 and woo = 4 x 0.3370, see Ref. 
These weights include the neutron scattering lengths of 
silicon and oxygen atoms; the factors of 2 and 4 account 
for the chemical composition of silica (two oxygen atoms 
for every silicon atom). 

To obtain the Fourier-broadened correlation function 
Tb(r), the correlation function of Eq. (||) is transformed 
to obtain the interference function Qi{Q): 



Qi{Q) - 



[T(r) — r] sm{Qr)dr, 



(6) 



r=0 



where Vmax is half the edge of the cubic simulation vol- 
ume of the model. Next, we transform the interference 
function back to obtain Ti,(r): 



Tb{r)=r+- 

TT 



M[Q)Qi{Q) sm{Qr)dQ, (7) 



using the experimentally relevant Qmax- The function 
M{Q) is an (optional) modification function commonly 
used to reduce Fourier artifacts. A popular choice is the 
Lorch modification function |18|] given by: 



M{Q) = 



Qr. 



ttQ 



ttQ 



(8) 



The neutron T{r) of vitreous silica was accurately mea- 
sured by Grimley, Wright and co-workers in 1991 [|l9[ |2^ 
using high-energy neutrons with Qmax = 45.2 . For 
each of the three thermalized networks we determined the 
Fourier- broadened correlation function Ti,{r) using the 
same value for Qmax and the Lorch modification func- 
tion. The experimental correlation function and the cor- 
relation function of the thermalized '60k' network are 
compared to each other in Fig. |[ We observe excellent 
agreement. As a more quantitative measure of the agree- 
ment between model and experiment we consider the 
factor |15|1 defined as: 



Tj2 / [Texpif) 



Tmodeiir)] dr 



(9) 



The discrepancy between the experimental T(r) and 
Tb{r) of the thermalized '60k' network, over the range 
1.0 < r < 10.0 A, equals = 5.1%. For the thermal- 
ized networks '3k' and '300k' we obtain R^ factors of 
4.7% and 4.9%, respectively, over the same range. 

In the top frame of Fig. ^ we compare the total in- 
terference function of the thermalized '60k' network to 
experiment. There is excellent agreement on the overall 
peak positions and also on the damping of the interfer- 
ence function for large Q. We observe some discrepancy 
at Q « 15 A~^, where the small peak visible in the ex- 
perimental data is not reproduced by our network. The 
discrepancy seems to be temperature-related: the lower 
frame of Fig. ^ shows that the network at zero tempera- 
ture does reproduce the peak correctly. 

We observe in Fig. ^ that the total correlation func- 
tion is significantly broadened by Fourier and thermal 
effects, in particular the first peak at r « 1.6 A. We have 
quantified the broadening of this peak by fitting it to the 
form: 



P{r) = 



1 



V 



'271 ra 



cxp 



{r - roY 
2ct2 



(10) 



with fit parameters {rp, cr, 77}. This follows the same pro- 
cedure used in Ref. |l^ to obtain bond lengths and coordi- 
nation numbers from experimentally obtained correlation 
functions. The parameters Tq and a provide estimates for 
the mean Si-O bond length and its variation. The pa- 
rameter r/ can be used to extract the mean coordination 
number z of silicon atoms: 

wsiSi + wsio + woo 



Airripo- 



WSiO 



(11) 



'=0 



with po the number density of oxygen atoms, Wap 
the neutron weights and T{r) normalized such that 
limr^oo T{r) = r. 
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FIG. 3; Total neutron correlation function T(r) of vitreous 
silica. The top frame shows the Fourier-broadened correla- 
tion function of thermalized network '60k' (solid) compared 
to the experimental result of Wright and co-workers (dashed) 
taken from Ref. ^ The lower frame shows the correlation 
function of the same network at zero temperature without 
Fourier broadening. 



FIG. 4: Total neutron interference function Qi{Q) of vitre- 
ous silica. The top frame shows Qi{Q) of thermalized net- 
work '60k' (solid) compared to experimental data (dashed) 
of Grimley and co-workers |l9j. The lower frame shows the 
interference function of network '6Qk' at zero temperature. 



The effects of Fourier and thermal broadening on the 
parameters {ro,(T, 2} are reported in Table III. As ex- 
pected, the broadening is significant: the combined ef- 
fect of Fourier and thermal broadening will boost cr by a 
factor of over five. 

It is also interesting to compare the silicon coordina- 
tion number obtained directly from the atomic coordi- 
nates of a network, to that obtained from a fit to the 
correlation function of that network. To this end the data 
of the networks at 300 K without Fourier broadening in 
Table [II are compared to Table ||. The fitting procedure 
accurately predicts the Si-0 bond length and its varia- 
tion (7, but not the coordination number which we know 
to be 4.000. This result is important because it shows 
that in the ideal situation free of Fourier broadening, fits 
based on Eq. (10) will systematically underestimate the 
coordination number. The cause of this underestimation 
is due to the non-gaussian nature of the peak in rT{r): 
fitting with a more elaborate function, for instance with 
a sum of two gaussians, yields the correct coordination 
number. The table also shows that fits with Eq. (|l^) ap- 
plied to Fourier-broadened data tend to overestimate the 
coordination numbers. 



TABLE III: The effect of Fourier broadening (FB) and ther- 
mal broadening on the optimal parameters {ro, cr, 2} obtained 
from a fit of Eq. ( |lo| ) to the first peak in T(r) of our networks. 
Network '3k' proved to be too small to accurately sample the 
first peak and was not used in this comparison. The param- 
eters ro and a are given in A. 







'60k' 


'300k' 


r = o K 


ro 


1.607 


1.601 


no FB 


a 


0.011 


0.011 




z 


3.962 ±0.011 


3.980 ± 0.004 


r = Q K 


ro 


1.609 


1.607 


with FB 


a 


0.048 


0.048 




z 


4.155 ±0.052 


4.161 ±0.053 


T = 300 K 


ro 


1.611 


1.610 


no FB 


a 


0.041 


0.040 




z 


3.977 ±0.013 


3.978 ±0.011 


T = 300 K 


ro 


1.614 


1.612 


with FB 


a 


0.060 


0.060 




z 


4.046 ±0.019 


4.048 ±0.019 
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VII. CONFIGURATIONAL ENTROPY 

Recently, we developed a method to determine the con- 
figurational entropy of a network pl} |. The method re- 
quires only the atomic coordinates and a list of bonds 
between particles, for a single well-relaxed configuration. 
In the original paper, the method was applied to a sil- 
ica network consisting of 3000 atoms. Starting with the 
atomic coordinates of this network, the list of bonds was 
constructed using rr = 1.80 A for the Si-0 bond-length 
cutoff. Each Si-O-Si bridge was replaced by a single Si-Si 
bond. We then applied our method to the resulting sil- 
icon backbone network to obtain for the configurational 
entropy s = 0.88 fcs per silicon atom. 

The accuracy of the method is sensitive to the size of 
the network. As mentioned in the original paper, the 
configurational entropy of vitreous silica reported above 
is most likely an underestimate due to the limited size 
of that network. The large networks generated in this 
work allow us to quantify these finite size effects. They 
also allow us to determine the minimum size of a network 
that is required to accurately measure the configurational 
entropy. 

To this end, we use our method to determine the en- 
tropy of networks '3k', '60k' and '300k' in two differ- 
ent ways: first, by using the entire simulation volume of 
the network; secondly, by using only half the volume of 
the simulation cell. The first procedure yields entropies 
of 0.83, 0.99 and 1.04 ks per silicon atom, respectively, 
while the second procedure yields 0.78, 0.97 and 1.03 fcs 
per silicon atom, respectively. This illustrates that finite- 
size effects are significant for the smallest network, but 



small for the large ones. Most of the difference in entropy 
between the three networks is due to the varying degree 
of relaxation in these networks. 



VIII. CONCLUSIONS 

In summary, we have presented a Monte Carlo-based 
approach for the generation of well-relaxed networks 
of vitreous silica, based on earlier work of Tu et al. 
With this method, networks containing 3000, 60,000 and 
300,000 atoms are generated. Compared to networks 
generated with molecular dynamics, our networks have 
smaller bond-angle variations and are nearly defect-free, 
indicating they are better relaxed. 

The total correlation function T(r) of our networks is 
in excellent agreement with neutron scattering experi- 
ments (with i?-^ factors around 5%), provided that ther- 
mal effects and the maximum experimental inverse wave- 
length Qmax are included in the comparison. We also de- 
termined the silicon coordination number of our networks 
by fitting a gaussian to rT(r), as is commonly done in ex- 
periments. We observe that the results of this procedure 
are biased to lower numbers by the non-gaussian nature 
of the peaks, and to higher numbers due to finite value of 
Qmax ■ The configurational entropy of vitreous silica was 
determined to be 0.99 ks and 1.04 fc^ per silicon atom, 
for the networks containing 60,000 and 300,000 atoms, 
respectively. 

Upon request, we will communicate the atomic coor- 
dinates of our networks. 
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